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PREFACE 
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WAVE  GROUP  ANALYSIS  BASED  ON  KIMURA'S  METHOD 


PART  I :  INTRODUCTION 

1.  High  sea  waves  tend  to  appear  in  groups  rather  than  individually. 
Engineers  are  finding  that  this  grouping  has  important  ramifications  on  the 
motions  and  resonances  of  moored  structures  and  vessels,  harbor  resonance, 
stability  and  overtopping  of  shore  protection  structures,  and  surf  beat. 
Because  of  the  nature  of  wave  grouping,  its  prediction,  control,  and  analysis 
are  especially  important  in  shallow-water  laboratory  basins  such  as  the 
Coastal  Engineering  Research  Center's  (CERC's)  directional  spectral  wave 
basin. 

2.  This  report  is  the  result  of  research  conducted  at  the  Hydraulics 
Laboratory  of  the  National  Research  Council  of  Canada  (NRCC) ,  Ottawa,  Ontario, 
Canada,  from  4-20  September  1985.  During  this  time,  the  original  version  of 
computer  program  KIMUR5  was  researched,  written,  debugged,  and  tested.  The 
program  calculates  wave  group  run  probabilities,  lengths,  means,  and  standard 
deviations  using  Kimura's  method  (Kimura  1980).  His  method,  which  is  based  on 
the  assumption  that  successive  wave  heights  are  mutually  correlated,  has  been 
demonstrated  (van  Vledder  1983b;  Thomas,  Baba,  and  Harish  1986)  to  be  superior 
to  Goda's  method. 

3.  This  report  describes  wave  grouping  and  the  differences  in  theory 
between  Goda's  and  Kimura's  methods.  Additionally,  it  documents  the  computer 
program  KIMUR5  and  serves  as  a  user's  manual  for  program  organization,  input/ 
output  operations,  and  test  cases.  Finally,  it  provides  recommendations  for 
future  expansion  of  the  program.  A  copy  of  the  computer  program  KIMUR5  and 
associated  subroutines  is  available  upon  request. 


PART  II:  WAVE  GROUPING 


4.  Bounded  long  waves  are  associated  with  the  occurrence  of  wave  groups 
and  produce  a  variation  of  the  mean  water  level  which  produces  a  setdown  under 
wave  groups  and  a  setup  between  groups  (Figure  1).  Longuet-Higgins  and 
Stewart  (1962)  first  described  this  second-order  or  nonlinear  effect  which 
results  from  a  variation  in  the  radiation  stress  (proportional  to  the  square 
of  the  local  wave  height).  The  forced  long  wave  propagates  at  the  group 
velocity  of  the  primary  waves.  Its  amplitude  is  proportional  to  the  square  of 
the  wave  envelope  and  is  relatively  small,  but  it  can  increase  dramatically  as 
the  depth  and  frequency  decrease  and  wave  groupiness  increases.  The  second- 
order  wave  system  propagates  with  phase  opposite  to  the  envelope  of  the  first- 
order  system.  A  crest  of  the  second-order  system  coincides  with  a  trough  of 
the  wave  group  envelope.  Second-order  currents  are  also  created  by  the  occur¬ 
rence  of  wave  groups.  These  currents  are  important  in  the  calculation  of 
resistance  forces  of  structures  and  mooring  forces  for  vessels. 

5.  A  succession  of  high  waves  that  exceeds  some  arbitrary  threshold 
value  (typically  median,  mean,  or  significant  wave  height)  is  called  a  run  of 
high  waves,  and  the  number  of  waves  in  this  run  is  the  run  length  (Figure  2). 
The  total  or  complete  run  is  the  combination  of  the  run  of  high  waves  followed 
by  the  run  of  iow  waves  (i.e.  succession  of  waves  which  fall  below  the  thresh¬ 
old  value).  The  total  run  is  analogous  to  the  zero-upcrossing  period  of  the 
wave  profile,  except  that  the  time  series  is  composed  of  individual  wave 
heights  rather  than  surface  elevations.  Reference  to  a  wave  group  assumes 
that  a  run  of  high  waves  is  intended. 


WAVE  ENVELOPE 
■  2ND-ORDE  R  WA  VE 
■  1ST-0RDER  OR  PRIMARY  WA  VE 
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Figure  1.  Schematic  of  wave  grouping  and  second-order  effects 


TOTAL  RUN  LENGTH 


NUMBER  OF  SUCCESSIVE  WAVES 

Figure  2.  Definition  of  wave  group  run  lengths 


PART  III:  THEORY  OF  WAVE  GROUPS 


Goda's  Method 


6 .  Goda ' s 

or  independent. 

probability  p* 

H  is 
c 


method  assumes  that  successive  wave  heights  are  uncorrelated 
The  derivation  is  based  on  probability  theory.  If  the 
that  a  wave  height  H  will  be  greater  than  a  threshold  value 


p  =  Prob  [H  >  H  ]  (1) 

c 

then  the  probability  that  H  will  be  less  than  or  equal  to  Hc  is  given  by 


q  =  Prob  [H£H]“l-p 


(2) 


since  p  +  q  =  1  .  (Some  authors  use  the  symbol  for  the  cutoff  or 

threshold  wave  height.) 

Run  of  high  wave  statistics 

7.  The  probability  distribution  of  a  run  of  successive  high  waves 

follows  the  geometric  distribution 


P(Jj) 


V1 

p  q 


j 1  *  1,2,3... 


(3) 


which  means  that  -  1  successive  waves  will  exceed  the  threshold,  and  the 

jj*1  wave  will  not.  The  probability  of  successive  wave  heights  exceeding  the 
threshold  is  equal  to  the  product  of  separate  probabilities  for  each  event 
because  the  wave  heights  are  independent  in  Goda's  (1985)  theory. 

8.  According  to  Goda  (1985),  the  mean  run  length  jj  and  standard 
deviation  of  the  length  of  a  run  of  high  waves  o(jj)  are,  respectively, 


*  For  convenience,  symbols  and  abbreviations  are  listed  in  the  Notation 
(Appendix  D) . 


00 


(4) 


h  ■  E(JlJ  ■  2 
V1 


q 


0(j,>  -  e[jJ]  -  E2[J,1  -  G- 


(5) 


where  E[  ]  is  the  expectation  operator. 

Total  run  statistics 

9.  For  a  total  run  of  successive  waves,  the  probability  distribu¬ 

tion  P ( j »  ®ean  ,  and  standard  deviation  oCjj)  are,  respectively, 

j_-l  Jo-* 

p ( j 2 )  “  ~  P  -  q  j2  =  2,3,4,...  (6) 


j2  =  £tj2^ 


o(j2)  *  E(j2] 


E2  [  j  2 1 


a_ 

2 

P 


(7) 

(8) 


Example  calculations 

10.  Example  calculations  based  on  Equations  1  to  8  for  the  probabilities 
and  mean  and  standard  deviations  for  the  run  of  high  waves  and  total  run  are 
listed  in  Table  l  for  various  values  of  threshold  wave  height.  According  to 
Coda's  model,  measured  average  group  length  larger  than  the  values  given  in 
Table  1  indicates  a  higher  level  or  degree  of  wave  group  formation. 


Kimura's  Method 


11.  Kimura's  model  assumes  that  successive  wave  heights  are  mutually 
correlated  and  form  a  Markov  Chain.  The  concept  of  mutual  correlation  implies 
that  successive  waves  are  dependent  or  correlated.  A  high  wave  rarely  appears 
by  itself;  rather,  it  is  more  likely  to  be  followed  by  other  high  waves.  It 
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Table  1 

Theoretical  Wave  Group  Statistics  for  Various  Wave  Height  Thresholds 

Goda's  Model 


Threshold 

Wave  Height 

Quantity 

Median 

Mean 

Significant 

Highest  1/10 

P 

0.500 

0.456 

0.135 

0.039 

q 

0.500 

0.544 

0.865 

0.961 

2.00 

1.84 

1.16 

1.04 

O(jj) 

1.41 

1.24 

0.42 

0.21 

h 

4.00 

4.03 

8.58 

26.55 

a(j2) 

2.00 

2.04 

6.92 

25.01 

seems  that  the  waves  have  a  "memory"  which  dictates  that  one  high  wave  will  be 
followed  by  another  high  wave  rather  than  a  low  wave. 

12.  Transition  probabilities  for  simultaneous  exceedance  and  non¬ 
exceedance  of  the  threshold  wave  height  are  calculated  based  on  the  ratio  of 
one-  and  two-dimensional  (i.e.  joint  or  bivariate)  Rayleigh  probability 
density  functions  (PDF).  From  these  Rayleigh-derived  transition  probabili¬ 
ties,  the  probability  of  a  run  of  various  lengths,  the  average  run  length,  and 
the  standard  deviation  of  the  run  length  are  calculated  for  a  run  of  succes¬ 
sive  high  waves  and  a  total  run. 

Markov  Chain 

13.  A  fundamental  assumption  of  Kimura's  model  is  that  successive  wave 
heights  form  the  Markov  Chain.  The  transition  equation  describing  the  Markov 
Chain  is 


P 

n 


V 


n 


where 

P  =  distribution  after  n-time  transitions 
n 

Pq  =  initial  distribution 
p  =  transition  probability  matrix 


(9) 


If  a  threshold  wave  height  Hc  Is  selected  (l.e.  mean,  median,  significant, 
or  highest  1/10  wave  height),  waves  with  height  H  will  fall  into  one  of  two 
states  or  groups  as  shown  below 

State  Condition 

1  H  <  H 

c 

2  H  >  H 

c 

The  initial  distribution  is  then 


PQ  =  (0,1) 


(10) 


since  a  run  of  high  waves  begins  when  State  2  is  first  reached.  The  transi¬ 
tion  probability  matrix  is  given  by 


P 


11 

p12 

21 

P22 

(ID 


where  the  individual  elements  or  conditional  probabilities  are  defined  as 


=  Prob  H. 

< 

H 

1  H- 

< 

H 

11 

i+1 

c 

1  i 

( 

=  Prob  IV  , 

> 

H 

I  H. 

< 

H 

12 

L  1+1 

c 

• 

=  Prob  k  , 

< 

H 

1  H. 

> 

H 

21 

L  1+1 

c 

1  i 

» 

=  Prob  k  , 

> 

H 

1  H . 

> 

H 

22 

L  i+1 

c  1 

1  i 

( 

(12) 


and  and  represent  successive  wave  heights.  Thus,  p^  is  the 

probability  that  neither  successive  wave  exceeds  the  threshold  Hc  ,  and  p ^ 
is  the  probability  of  simultaneous  exceedance  by  both  wave  heights.  By  sub¬ 
stituting  Equations  10  and  11  into  Equation  9,  processing  for  n-tirae  transi¬ 
tions  and  precluding  the  transition  probabilities  from  State  1,  we  obtain  the 
probability  distribution  for  the  run  of  high  waves.  If  simple  induction  con¬ 
tinues,  the  probability  distribution  of  the  total  run  can  be  determined  in  a 
similar  fashion. 
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Transition  probabilities 

14.  The  other  fundamental  assumption  of  Kimura's  method  is  that  the 
transitional  or  conditional  probabilities  for  wave  height,  p^  and  p 22  , 
are  defined  in  terms  of  the  Rayleigh  distribution.  The  conditional  probabili¬ 
ties  are 


H  H 
c  c 


/  f  p(H1,H2)dH1< 


0  0 

*n  * - if 


/  q(H1)dH1 


Jj 

+*■' 

‘ft 

u 

i 

m 

•f 


00  00 


f  f  p(H1,H2)dH1< 


H  H 
c  c 


/  q(H1)dH1 


jft 

to 


where 

p(Hj,H2)  ■  joint  or  two-dimensional  Rayleigh  PDF  for  two  successive 
wave  heights 

Hi,H2  *  dummy  wave  height  variables 

q(Hj)  *  Rayleigh  PDF  for  individual  wave  heights 

15.  The  PDF  q(Hj)  defined  in  terms  of  the  root-mean-square  wave  height 

H  and  the  mean  wave  height  H  is 
r  in 


q  (H,  )  -  -a-  exp  I - - 


H, 

ir  1 

~  — J  exP 

‘  H2 
m 


•‘V* 


16.  The  Rayleigh  joint  PDF  is  similarly  defined  by  Kimura  (1980)  based 
on  earlier  work  of  Rice  (1944,  1945)  and  Uhlenbeck  (1943)  as 


.< . 

» <  *| 


Ld 


u* 


TT. 


V  !*  I  . 


*»  L 


»» 

/i 

> 


1 

ji 


p(HrH2) 


4HjH2 

(i  -  <2)h* 


,2H1H2 
4(1  -  *c2)H^ 


2  2 
Hj  +  H2 


(1  -  O 


2  2 
(1  -  <) 


2  2 
H1  +H2 


4(1  -  O 


2  2 
2d  -  <)  ir 

m 


where  <  Is  the  correlation  parameter  and  Iq[  ]  is  a  modified  Bessel  func¬ 
tion  of  zeroth  order. 

Correlation  parameter 

17.  To  solve  for  the  joint  Rayleigh  PDF  and  the  associated  transition 
probabilities  j  and  p22  ,  the  correlation  parameter  <  is  required. 
Kimura  (1980)  and  some  authors  define  it  in  terms  of  the  variable  P  as 


<  *  2p 


Uhlenbeck  (1943)  showed  that  the  correlation  parameter  <  is  related  to  the 
correlation  coefficient  R^O)  ,  a  measure  of  the  degree  of  correlation  or 
dependence  between  successive  wave  heights,  by 


Rhh(1) 


(1  -  O  K(k) 
_ 2 _ 


where  E(  )  and  K(  )  are  complete  elliptic  integrals  of  the  first  and 
second  kind,  respectively.  Some  authors  use  Y^  to  define  the  correlation 
coefficient.  The  correlation  coefficient  range  is  0  to  1.0.  A  value  of  zero 
corresponds  to  the  Goda  model.  Several  investigators  (Goda  1985)  have  calcu¬ 
lated  correlation  coefficients  of  0.24  for  successive  wind  waves  and  0.5  -  0.8 
for  swell.  The  amount  of  correlation  tends  to  increase  with  higher  wave 
heights  and  narrower  wave  spectra. 

18.  Battjes  (1974)  demonstrated  that  an  infinite  series  representation 
for  the  elliptic  integrals  could  be  used  to  approximate  Equation  18  as 


Rhh(1)  =  4(4  -  it) 


(19) 


If  this  relation  is  inverted,  the  correlation  parameter  <  can  be  determined 
given  the  correlation  coefficient  R,,  (1)  as  follows: 


i 

!  where 


(20) 


R  = 


(21) 


Battjes  (1974)  found  that  this  approximation  is  very  good  for  correlation 
coefficients  less  than  0.7  to  0.8.  From  0.8  to  1.0  the  difference,  although 
slight,  is  still  noticeable. 


Methods  for  deter- 

mining  correlation  parameter 

19.  The  correlation  parameter  can  be  determined  in  one  of  four  ways: 

a.  Time  Domain  Method  1:  assumed  correlation  coefficient. 

b.  Time  Domain  Method  2:  autocorrelation  technique. 

£.  Frequency  Domain  Method  1:  Goda's  spectral  peakedness 
parameter. 

d.  Frequency  Domain  Method  2:  Battjes'  spectral  derivation. 

20.  Assumed  correlation  coefficient.  The  assumed  correlation  coeffi¬ 
cient  method  is  the  one  currently  coded  in  the  computer  program  KIMUR5.  The 
input  value  of  the  correlation  coefficient  R^^(l)  is  based  on  field  measure¬ 
ments  for  similar  wave  conditions  as  the  wave  height  time  series  to  be 
analyzed.  The  correlation  parameter  <  is  determined  indirectly  using  Equa¬ 
tions  19  and  20  above. 

21.  Autocorrelation  technique.  In  the  autocorrelation  technique  the 
correlation  coefficient  R^(k)  is  first  calculated  from  a  zero-meaned,  mea¬ 
sured,  or  simulated  wave  height  time  series.  The  correlation  parameter  ic  is 
determined  indirectly  using  Equations  19  and  20  as  before.  The  autocorrela¬ 
tion  function  estimate  is  normalized  by  the  variance  of  the  wave  height  time 


series  to  give  the  correlation  coefficient  defined  as 


Rhh(k) 


J.  1 

2  N  -  k 

o.. 


N-k 

I 


i=l 


H.HU. 

1  i+k 


k  *  1,2,3,... 


(21) 


where  N  is  the  total  number  of  points  in  the  wave  height  time  series,  and 
o  is  the  standard  deviation  of  the  series.  The  lag  k  is  the  difference  in 
number  between  wave  heights  and  is  equal  to  1  for  successive  wave  heights. 

For  every  other  wave  height,  the  correlation  coefficient  would  be  written  as 
R^h(2)  ,  every  third  wave  height  R^(3)  »  etc.  The  dependency  between  wave 
heights  has  been  found  by  several  investigators  (van  Vledder  1983a)  to 
decrease  rapidly  as  the  lag  is  increased  beyond  successive  wave  heights  (i.e. 

Rhh(1))* 

22.  Goda's  spectral  peakedness  model.  The  spectral  peakedness  model  is 
a  frequency  domain  model  based  on  a  relationship  between  wave  grouping  and 
spectral  form  investigated  by  Goda  (1970,  1976),  Yamaguchi  (1981),  and  Kimura 
(1980)  among  others.  It  is  based  on  the  spectral  peakedness  factor 
defined  by  Goda  as 


Q 


P 


f  S2(f)  df 


(22) 


where 

m  =  zeroth  moment  of  the  time  series 

o 

f  ■  frequency 

S(f)  =  spectral  estimate  of  the  surface  elevation 
Goda  (1985)  found  the  spectral  peakedness  parameter  to  be  insensitive  to  the 
high  frequency  cutoff  used  in  spectral  analysis.  Its  value  ranges  between  1 
for  white  noise,  2  for  wind  waves,  and  U  to  8  for  swell  conditions.  The 
investigators  mentioned  above  found  that  the  average  group  length  increases  as 
the  peakedness  parameter  increases,  and  a  narrow  spectrum  has  a  greater  degree 
of  grouping  than  a  widebanded  spectrum. 

23.  Based  on  field  wave  data  and  numerical  simulations  for  large  values 


of  Q  ,  Ewing  (1973)  proposed  an  approximately  linear  relationship  between 
the  mean  run  length  and  the  spectral  peakedness  parameter  for  a  given 

cutoff  or  threshold  wave  height  Hc  as  follows: 


r-i 

J1  H 

c 


'2m 


(23) 


24.  Goda  (1970)  proposed  the  relationship  between  the  correlation 
and  the  wave  peakedness  parameter  Q  as  shown  in 


coefficient 


Rhh(1) 


Figure  3  (obtained  from  numerical  simulations) . 


Figure  3.  Relation  between  correlation  coefficient  and 
spectral  peakedness  parameter 


25.  Battjes*  spectral  derivation.  Based  on  earlier  work  of  Arhan  and 
Ezraty  (1978)  on  correlations  with  joint  PDF’s  and  Rice  (1944,  1945)  on 
theoretical  envelope  statistics,  Battjes  and  van  Vledder  (1984)  showed  that 
the  correlation  parameter  <  can  be  calculated  spectrally  by 
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[  S(f)  cos  (2irf  T  )  df 
J  m 
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[  S(f)  sin  C2irf  T  )  df 
J  m 


(24) 


where  is  the  mean  wave  period  obtained  from  zero-crossing  or  spectral 

analysis.  In  this  sense,  the  correlation  parameter  <  is  a  measure  of  the 
spectral  width,  and  Battjes  and  van  Vledder  (1984)  noted  that  it  is  more 
"robust"  than  Goda's  spectral  peakedness  parameter  since  it  is  not  biased 

by  sampling  variability. 

Run  of  high  wave  statistics 

26.  In  Kimura's  model,  the  probability  of  a  run  of  successive  high 
waves  of  run  length  j  ^  is  defined  in  terms  of  the  transition  probability 
P22  as 


V1 

p (J  1 )  =  P22  (1  -  P22) 


1,2,3, 


(25) 


27.  The  mean  jj  and  standard  deviation  o(jj)  are,  respectively. 


and 


J>  ‘  0  -‘p22> 


o(jj) 


VP  22 

(1  -  P22) 


(26) 


(27) 


Similarity  exists  among  Equations  25  to  27  and  Equations  3  to  5  for  Goda's 
method  in  which  p^^  and  (1  -  p^2)  replace  p  and  q  ,  respectively. 

Total  run  statistics 

28.  For  a  total  run,  the  probability  distribution  P(j2)  ,  mean  , 
and  standard  deviation  o(j2)  are,  respectively, 


P(j2) 


(1  -  Pu)(l  -  P22) 


j2_1 
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12  -  2,3,4,...  (28) 
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(29) 
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[(1  -  p22)^  (1  -  pn)ZJ 


(30) 


Again,  there  is  similarity  with  the  total  run  statistics  defined  for  Goda's 
model  in  Equations  6  to  8. 


Comparison  of  Methods 

29.  Goda's  model  for  wave  group  run  lengths  assumes  that  wave  heights 
are  independent  or  uncorrelated,  although  Rayleigh  distributed.  Rye  (1974), 
Kirnura  (1980),  and  others  have  shown  that  wave  heights  are  positively  corre¬ 
lated.  Thus,  in  comparisons  with  actual  field  measurements  for  varying  wave 
environments  (including  wind  wave  generation  in  storms)  by  several  investiga¬ 
tors,  Goda's  method  yields  a  constant  value  for  several  values  of  the  corre¬ 
lation  coefficient  that  seriously  underpredict  the  degree  of  wave  groupings. 
These  comparisons  of  field  measurements  with  Goda's  model  values  for  median 
Hmed  an<*  si£rif*cant  Hs  wave  height  threshold  values  are  listed  in  Table  2 
(van  Vledder  1983a)  which  shows  that  the  measured  values  for  the  run  lengths 
are  greater  than  those  Coda  predicted. 

30.  Table  3  shows  the  results  of  some  computer  simulations  by  Kirnura 
(1980)  for  the  group  lengths  of  a  run  of  high  waves  for  spectra  of  various 
peakedness  and  uniform  phase  distributions.  Goda's  and  Kimura’s  theoretical 
values  for  five  different  correlation  coefficients  are  compared  with  the 
simulated  data  for  threshold  wave  heights  equal  to  the  mean  and  significant 
wave  height.  Kimura's  model  shows  a  strong  agreement  with  the  data,  while 
Goda's  model  gives  a  constant  value  that  underpredicts  the  group  length. 

31.  Table  4  contains  analogous  results  by  Goda  (1983)  for  the  group 
lengths  of  a  run  of  high  waves  using  measured  data  representative  of  long 
traveled  swell  with  a  narrow  spectrum  and  high  correlation  coefficients. 
Again,  there  is  serious  underpredict  ion  of  the  Goda  model  and  the  reasonable 
cor respondence  between  actual  and  predicted  run  lengths  using  the  Kirnura 
mode  1  . 


Table  2 

Comparison  of  Measured  Average  Group  Lengths  With  Goda’s  Model 

Run  of  High  Waves 
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Table  3 

Comparison  of  Goda's  and  Kimura's  Models  with  Simulated  Data 


Average  Group  Lengths  for  Run  of  High  Waves 


0.19 

0.23 

0.29 

0.33 

0.38 


Mean 

Klmura 


Threshold  Wave  Height 
Simulated  Goda 


2.20  1.15 
2.29  1.15 
2 . 34  1.15 
2.42  1.15 
2.45  1.15 
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Table  4 

Comparison  of  Goda*s  and  Kimura's  Models  with  Measured  Data 
Average  Group  Lengths  for  Run  of  High  Waves 


0.630  1.84  3.50  3.77  1.15  2.08  2.02 
0.688  1.84  3.84  4.15  1.15  2.29  2.49 
0.694  1.84  3.89  4.42  1.15  2.31  2.21 
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PART  IV:  COMPUTER  PROGRAM  DESCRIPTION 


32.  This  section  presents  the  documentation  for  the  main  program  KIMUR5 
and  the  eleven  associated  subroutines  required.  The  intent  is  to  provide  the 
user  with  the  documentation  necessary  to  run  the  program.  Appendix  A  contains 
a  listing  of  the  program  and  subroutines.  Appendix  B  contains  a  listing  of 
the  symbols  used  in  the  computer  program.  Appendix  C  lists  definitions  of 
parameters  for  each  subroutine  and  gives  other  documentation  including 
descriptions,  calling  statement,  calls  to  and  by  the  subroutine,  and 
references. 


Program  Specifications 


33.  Table  5  summarizes  the  program  specifications  for  program  KIMUR5. 

Table  5 

Program  Specifications  for  Program  KIMUR5 


Specification 

Description 

Computer 

DEC  VAX  11/750 

Location 

USAE  Waterways  Experiment  Station,  CERC 

Operating  system 

VAX/VMS  version  4.2 

Language 

Fortran  77 

Structure 

Interactive,  modular,  top  down 

Documentation 

Self-document ing 

Subroutines 

Eleven,  shelf-contained 

(Subroutines  BESI  &  QSF  from  NRCC  Scientific 

Library) 

Input 

Interactive  with  prompts,  logical  unit  5 

Output 

Disk  File  KIMUR.OUT,  logical  unit  2 

Accuracy 

Single  precision 

(Subroutine  BESI  requires  double  precision) 

Operating  procedure 

Compile:  FORTRAN  KIMUR5 

Link:  LINK  KIMVR5 

Run:  RUN  KIMUR5 

Input:  Enter  5  Input  values  at  keyboard 

Output:  TYPE  or  PRINT  KIMUR.OUT 

19 


Solution  Procedure 


34.  The  computer  code  is  presently  in  the  form  of  a  main  driver  program 
and  eleven  associated  subroutines.  Figure  4  is  a  flowchart  illustrating  the 
basic  steps  involved  in  program  calculations.  Table  6  lists  the  hierarchy  of 
the  individual  subroutines  in  the  program.  A  brief  description  of  each 
subroutine  is  contained  in  Table  7.  The  steps  indicated  in  Figure  4  ( 1 — D 
=  one  dimensional;  2-D  =  two  dimensional)  are  described  in  the  paragraphs 
below. 
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Input  Parameters 
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Correlation  Parameter 
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RPDF 
1-D  Rayl 
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Output  Result 
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Figure  4.  Flowchart  of  Program  KIMUR5 


20 


,  .  •/  s«.V  * '  VjV. 


r 

Table  6 

Hlerachy  of  Program  KIMUR5 


Main  _ Subroutines 

Program  Level  1  Level  2 

KIMUR5  INPUT 

KAPPA 

RPDF 

RJPDF  BESI 

RTP  RTPI 

HRUN 

TRUN 

OUTPT 


Level  3 


QSF 


Table  7 

Description  of  Subroutines  In  Program  KIMUR5 

Name  _ Description _ 

INPUT  Queries  user  for  input  parameters 

KAPPA  Calculates  correlation  parameter  K  given  correlation  coefficient 
RHH(l)  using  series  approximation  method  of  Battjes 

RPDF  Calculates  1-D  Rayleigh  PDF  Q(Hi)  for  individual  wave  heights 

RJPDF  Calculates  2-D  Rayleigh  joint  PDF  P(H1,H2) 

BESI  Calculates  Bessel  function  I  of  zeroth  order  (NRCC  Scientific  Sub¬ 
routine  Library) 

RTP  Calculates  Rayleigh  transition  probabilities  Pll  and  P22 

RTPI  Integrates  1-D  and  2-D  Rayleigh  PDF's  Q(H1)  and  P(H1,H2) 

QSF  Computes  vector  integral  values  for  a  given  equidistant  table  of  func¬ 

tion  values  using  combination  of  Simpson's  and  Newton's  3/8  Rules 
(NRCC  Scientific  Subroutine  Library) 

HRUN  Calculates  run  of  high  wave  group  statistics  of  probability  of  dif¬ 
ferent  run  lengths,  mean  run  length,  and  standard  deviation  of  run 
length 

TRUN  Calculates  total  run  group  statistics  of  probability  of  various  run 
lengths,  mean  run  length,  and  standard  deviation  of  run  length 

Outpt  Outputs  results  to  disk  file  for  display  and  archival 


Correlation  parameter 

35.  The  first  step  in  the  solution  procedure  is  the  calculation  of  the 
correlation  parameter  K  from  the  input  correlation  coefficient  RHH1.  Sub¬ 
routine  KAPPA  performs  this  operation  using  the  infinite  series  approximation 
for  the  elliptic  integrals  given  in  Equation  20. 

Rayleigh  probability  density  function 

36.  The  second  step  is  the  calculation  of  the  1-D  Rayleigh  PDF  by  Sub¬ 
routine  RPDF.  The  Rayleigh  PDF  given  in  Equation  15  for  the  mean  wave  height 
is  programmed  with  the  dummy  wave  height  variable  HI  defined  as 


HI 


N  *  DELH 


N  =  0,1,2, ...NH 


(31) 


Descriptions  of  the  symbols  used  in  the  computer  program  are  contained  in 
Appendix  B. 

Joint  Rayleigh 

probability  density  function 

37.  The  next  step  in  the  solution  procedure  is  the  calculation  of  the 
joint  Rayleigh  PDF  by  Subroutine  RJPDF  using  the  mean  wave  height  form  of 
Equation  16.  For  ease  of  programming,  it  is  calculated  in  terms  of  three 
factors  as 


P(H1,H2)  -  A  *  B  *  C 


(32) 


where  the  A  factor  is  a  constant  term,  the  B  factor  is  an  exponential  term, 
and  the  C  factor  is  the  modified  Bessel  function  of  zeroth  order.  Again, 
dummy  wave  height  interval  variables  HI  and  H2  are  used  and  defined  in  the 
range 


HI  =  N  *  DELH 


H2  =  N  *  DELH 


0,1,2, ...NH 


(33) 


38.  The  modified  Bessel  function  of  zeroth  order  is  evaluated  in  Sub¬ 
routine  BESI,  which  was  obtained  from  the  NRCC  Scientific  Subroutine  Library. 
It  was  verified  on  several  test  cases  using  a  Chemical  Rubber  Company  Handbook 
of  Mathematical  Sciences  (1978). 


Rayleigh  transition  probabilities 

39.  The  fourth  step  is  the  evaluation  of  the  Rayleigh  transition 
probabilities  Pll  and  P22  using  Equations  13  and  14,  respectively.  Subroutine 
RTP  sets  up  the  proper  lower  and  upper  array  element  integration  limits  of  the 
1-D  and  2-D  Rayleigh  PDF's  (i.e.  NL  and  NU,  respectively)  for  the  particular 
cutoff  (i.e.  threshold)  wave  height  HC  selected.  For  the  Pll  transition 
probability,  the  lower  and  upper  array  elements  are,  respectively, 


NL  =  1 


NU 


HC 

DELH 


(34) 


Similarly,  for  the  P22  transition  probability,  these  limits  are 


NL 


HC 

DELH 


NU  =  NH  +  1 


(35) 


40.  Subroutine  RTPI  is  called  by  subroutine  RTP  and  calculates  either 
transition  probability  Pll  or  P22  given  the  proper  lower  NL  and  upper  NU  array 
element  integration  limit.  The  discrete  1-D  and  2-D  Rayleigh  PDF's  are  inte¬ 
grated  numerically  using  Subroutine  QSF,  obtained  from  the  NRC  Scientific  Sub¬ 
routine  Library.  A  dummy  1-D  array  at  equidistant  points  is  evaluated  using  a 
combination  of  Simpson's  and  Newton's  3/8  rules.  It  has  been  thoroughly 
tested  by  NRC. 

Statistics  of  run  of  high  waves 

41.  The  fifth  step  is  the  calculation  of  group  statistics  for  a  run  of 
high  waves  based  on  the  transition  probability  P22.  The  probabilities  for 
various  run  lengths  PHR  are  calculated  using  Equation  25.  The  mean  run  length 
JIM  and  standard  deviation  of  run  length  SIGJ1  are  given  by  Equations  26 

and  27. 

Statistics  of  total  run 


42.  The  final  step  is  the  calculation  of  group  statistics  for  a  total 
run  using  transition  probabilities  Pll  and  P22.  The  probability  distribution 


PTR,  the  mean  run  length  J2M,  and  the  standard  deviation  of  the  run  length 
SIGJ2  are  defined  in  Equations  28,  29,  and  30,  respectively. 


Input  and  Output  Variables 


Input  variables 

43.  Subroutine  INPUT  queries  the  user  for  the  five  variables  listed  in 
Table  8.  Figure  5  is  an  example  of  the  input  required.  The  value  of  RHH1  is 
dimensionless  and  should  be  between  zero  and  unity.  (See  Part  III  of  this 
report  for  range  of  values  for  typical  wave  conditions.)  If  an  actual  time 
series  of  wave  height  measurements  is  used,  then  the  parameter  NH  should  be 
equal  to  the  total  number  of  points  in  the  series.  Otherwise,  a  value  of  NH 

Table  8 

Input  Variables 

_ Description _ 

Correlation  coefficient 

Total  number  of  wave  height  measurements  or  total  number  of  inter¬ 
vals  of  dummy  wave  height  variable  HI  and/or  H2 

Wave  height  increment  between  successive  HI  or  H2  wave  heights 

Mean  wave  height 

Threshold  wave  height 


Variable 
RHH 1 
NH 

DELH 

HM 

HC 


*  RUN  K I  MUR 5 

ENTER  1  FOR  TEST  CASE: 


ENTER  UAL UES  FOR: 

RHH1  CORRELATION  COEFFICIENT 
NH  TOTAL  #  OF  NAVE  HEIGHT  MEASUREMENTS 

DELH  WAVE  HEIGHT  INCREMENT 

UPPER  WAVE  HGHT  INTEGRATION  LIMIT  =  NH  *  DELH 
HM  MEAN  WAVE  HEIGHT 

HC  THRESHOLD  WAVE  HEIGHT 

.68  400  .0025  .33  .33 

FORTRAN  STOP 

Figure  5.  Example  input  format  for  Program  KIMUR5 


of  400  has  been  found  Co  give  reasonable  results.  The  program  is  dimensioned 
for  up  to  500,  however.  The  third  input  variable,  DELH,  corresponds  to  the 
width  of  the  class  interval  in  a  nondimensionalized  histogram  or  distribution 
function  of  wave  heights.  The  smaller  this  value,  the  more  accurate  the 
results.  The  product  of  NH  and  DELH  gives  the  upper  wave  height  integration 
limit  HU  which  should  be  greater  than  or  equal  to  the  largest  wave  height  in 
the  time  series.  The  mean  wave  height  should  be  determined  from  the  time 
series  if  actual  wave  heights  are  used.  According  to  Goda  (1985),  the 
relationship  between  the  maximum  wave  height  H  ^  and  the  mean  wave  height 
is  approximately 

H  =  0.31  to  0.39  H  (36) 

m  max 

since  the  significant  wave  height  H  =  1.6  H  and  H  *  1.6  to  2.0  H  , 

s  m  max  s 

depending  on  the  number  of  points  in  the  time  series.  A  value  for  H  of 

m 

0.33  is  representative  of  a  Rayleigh  distributed  wave  height  for  H  =  1.0  , 

max 

assuming  the  statistically  derived  maximum  wave  height  is  equal  to  the  largest 
wave  in  the  time  series  of  wave  heights.  Finally,  the  threshold  wave  height 
can  be  equal  to  the  mean,  median,  or  significant  wave  height.  Formulas 
relating  these  three  parameters  are 


H  «  0.939  H 
med  m 


H  =  1.597  H 
s  m 


(37) 


Output  variables 

44.  Output  variables  are  written  by  subroutine  OUTPT  to  a  disk  file  for 
later  viewing  or  printing.  Figure  6  is  an  example  of  the  output  file 
KIMUR.OUT.  The  five  input  variables  are  listed  along  with  the  upper  wave 
height  integration  limit.  The  correlation  parameter  K  and  the  two  transition 
probabilities  Pll  and  P22  are  written  in  the  "output  variables"  section.  For 
a  run  of  high  waves,  the  probabilities  PHR  of  a  run  of  successive  high  waves 
of  run  lengths  of  1  through  25  (in  10F7.3  format),  the  mean  group  length,  JIM, 
and  the  standard  deviation  of  the  group  length,  SIGJ1,  are  given.  The  first 
element  of  PHR  corresponds  to  the  probability  of  a  run  of  length  1,  the  second 
element  is  the  probability  of  a  run  of  2  waves,  the  third  a  run  of 


RESULTS  FROM  PROGRAM  KIMUR5 
GROUP  RUN  LENGTH  STATISTICS 


★★★★★INPUT  VARIABLES***** 


CORRELATION  COEFFICIENT,  RHH1  =  0.6800 
TOTAL  ♦  OF  WAVE  HEIGHT  MEASUREMENTS,  NH  *  400 
WAVE  HEIGHT  INCREMENT,  DELH  =  0.0025 
UPPER  WAVE  HGHT  INTEGRATION  LIMIT  =  1.0000 
MEAN  WAVE  HEIGHT,  HM  *  0.3300 
THRESHOLD  WAVE  HEIGHT,  HC  *  0.3300 


★★★★★OUTPUT  VARIABLES***** 


CORRELATION  PARAMETER,  K  = 

PROBABILITY  NEITHER  HI  NOR  H2  EXCEEDS,  Pll 
PROBABILITY  BOTH  HI  &  H2  EXCEED,  P22  = 


0.8399 
0 . 761? 
0.7202 


★★★★★HIGH  WAVE  RUN  GROUP  RESULTS***** 
PROBABILITIES,  PHR : 


0.280  0.202  0.145  0.105 

0.011  0.008  0.005  0.004 

0.000  0.000  0.000  0.000 

0.075 
0.003 
0 .000 

0 .0  54 
0.002 

0.039  0.028 

0.001  0.001 

0 .020 

0 .001 

0 .015 
0.001 

MEAN  GROUP  LENGTH,  JIM  = 

STD  DEV  GROUP  LENGTH,  SIGJ1 

= 

3.5743 

3 . 0334 

★★★★★TOTAL  WAVE  RUN  GROUP  RESULTS***** 
PROBABILITIES,  PTR : 

0.000  0.067  0.099  0.110  0.109  0.101 
0.045  0.037  0.030  0.024  0.019  0.015 
0.005  0.004  0 . 003  0 . 002  0 . 002 

0 .0  90  0.0  78 

0.012  0.010 

0 . 0  66 

0 . 0  0  8 

0 .0  55 
0 . 0  0  6 

MEAN  GROUP  LENGTH,  J2M  = 

STD  DEV  GROUP  LENGTH,  SIGJ2 

= 

7.7715 

4.7561 

Figure  6.  Example 

output 

format  from  Program  KIMUR5 

3  waves,  etc.  Similarly,  for  the  total 

run,  the 

first  25  probabilities 

PTR  of 

a  run  of  length  I  through  25,  the  mean  length,  J2M,  and  the  standard  deviation 
of  the  run  length,  SIGJ2,  are  given. 
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PART  V:  PROGRAM  VERIFICATION 


45.  In  this  section,  verification  of  the  program  with  sample  test  cases 
is  described  and  discussed. 


Test  Case  1 


46.  The  first  test  case  is  based  on  actual  wave  height  data.  Goda 
(1985)  describes  an  example  of  97  waves  with  a  mean  wave  height  of  2.1  m.  The 
wave  heights  are  distributed  in  a  range  from  0.1  to  5.5  m.  Thus,  the  inputs 
to  Program  KIMUR5  are  RHH 1  -  0.68,  NH  *■  97,  DELH  *  0.0567  (HU  *  NH  *  DELH 
-  5.5  m) ,  and  HM  «  2.1.  Table  9  summarizes  the  results  for  three  test  cases 
for  each  of  the  mean,  median,  and  significant  wave  heights  as  threshold  wave 
height  (using  Equation  37).  The  effect  of  various  threshold  wave  heights  on 
the  value  of  the  run  length  calculated  is  readily  apparent.  The  calculated 
values  appear  to  be  reasonable  when  compared  with  Goda's. 

Table  9 

Summary  of  Test  Case  1  Results 


Threshold  Wave 

Height 

Transition 

Probabilities 

Mean 

Run  Lengths 

Description 

HC 

Pll 

P22 

High 

Total 

Median 

2.0 

0.727 

0.735 

3.775 

7.431 

Mean 

2.1 

0.753 

0.719 

3.554 

7.594 

Signif icant 

3.3 

0.922 

0.543 

2.186 

14.922 

Test  Case  2 

47.  Van  Vledder  (1983b)  calculated  the  transition  probabilities  PI  1  and 
P22  and  the  mean  run  lengths  JIM  and  J2M  for  various  values  of  the  correlation 
coefficient  RHH1  and  threshold  wave  heights  of  the  median,  mean,  and  signifi¬ 
cant  wave  heights.  For  a  nondimensional  maximum  wave  height  of  1.0,  the  fol¬ 
lowing  inputs  were  used  in  program  KIMUR5:  RH1  «  .68,  NH  *  400,  DELH 
=  0.0025,  and  HM  =  0.33.  Table  10  shows  the  comparison  between  the  KTMUR5 
calculated  values  and  the  van  Vledder  values  (given  by  van  Vledder  to  four 
significant  places).  The  average  percent  error  listed  is  the  average  of  the 


Table  10 

Summary  of  Test  Case  2  Results 


Threshold  Wave  Height 

Transition 

Probability 

PI1  P22 

Mean  Run 
JIM 

Length 

J2M 

Average 

Percent 

Error 

Descript  ion 

Source 

H 

c 

Median 

van  Vledder 
KIMUR5 

0.31 

0.7371 

0.7334 

0.7371 

0.7378 

3.8037 

3.8146 

7.6073 

7.5650 

0.08 

Mean 

van  Vledder 
KIMUR5 

0.33 

0.7651 

0.7617 

0.7197 

0.7202 

3.5674 

3.5743 

7.8243 

7.7715 

0.35 

Significant 

van  Vledder 
KIMUR5 

0.53 

0.9314 

0.9315 

0.5592 

0.5525 

2.2686 

2.2348 

16.8359 

16.8360 

0.67 

absolute  values  of  the  four  errors  between  the  two  transition  probabilities 
and  the  two  mean  run  lengths.  These  errors  are  defined  as  the  difference 
between  van  Vledder's  and  the  calculated  values  for  each  quantity  divided  by 
van  Vledder's  value.  Thus,  the  average  percent  error  is  less  than  0.67  per¬ 
cent  for  all  cases  tested  and  shows  very  good  agreement  with  van  Vledder's 
results . 

48.  Van  Vledder*  recommends  that  the  total  number  of  increments  NH 
should  be  greater  than  150  to  200.  Table  11  lists  the  differences  in 
calculated  values  for  a  threshold  wave  height  equal  to  the  mean  wave  height 
for  various  NH  values  of  50,  100,  200,  400,  and  500.  The  average  percent 
error  between  the  program's  values  and  van  Vledder's  decreases  markedly  for 
increases  in  the  number  of  intervals. 

Discussion  of  Results 

49.  The  agreement  of  the  KIMUR5  model  with  van  Vledder's  calculated 
values  is  excellent.  Many  factors  could  account  for  the  slight  differences 
observed.  According  to  van  Vledder,*  his  program  is  dimensioned  for  double 
precision,  explicitly  calculates  the  1-D  Rayleigh  integral,  and  uses  a 


*  Personal  Communication,  2  January  1986,  with  Dr.  G.  Ph.  van  Vledder,  Delft 
University  of  Technology,  Department  of  Civil  Engineering,  Delft,  The 
Netherlands . 


Table  11 

Effect  of  Various  NH  Values 


NH 

Transition 

Probabilities 

Mean  Run  Lengths 

Average 

Percent 

Error 

Pll 

P22 

High 

Total 

van  Vledder 

0.7651 

0.7197 

3.5674 

7.8243 

— 

50 

0.7219 

0.7445 

3.9141 

7.5105 

5.70 

100 

0.7514 

0.7268 

3.6605 

7.6837 

1.79 

200 

0.7584 

0.7224 

3.6025 

7.7408 

0.83 

400 

0.7617 

0.7202 

3.5743 

7.7715 

0.35 

500 

0.7624 

0.7198 

3.5687 

7.7778 

0.25 

maximum  wave  height  value  of  10  times  the  mean  wave  height  as  the  upper  wave 
height  integration  limit.  His  program  also  iteratively  checks  for  the  optimum 
number  of  steps  to  use  in  calculating  the  values  of  the  joint  Rayleigh 
integral . 

50.  My  investigations  showed  that  double  precision  did  not  make  any 
noticeable  difference  (to  E-04  precision  for  the  output  results)  for  a  value 
of  NH  of  100.  Different  programming  techniques,  the  numerical  integration 
routine  used,  explicit  calculation  of  the  1-D  Rayleigh  integral,  and  the 
method  of  selection  of  the  lower  and  upper  cutoff  limits  in  the  integrals 
probably  account  for  the  slight  differences  observed. 


PART  VI:  SUMMARY  AND  RECOMMENDATIONS 


51.  The  coastal  engineering  research  community  has  recognized  the  need 
to  model  wave  groups  as  well  as  spectral  waves.  High  waves  in  groups  can  pro¬ 
duce  more  damage  than  isolated  high  waves,  and  engineers  are  finding  that  this 
groupiness  has  important  ramifications  in  the  motions  and  resonances  of  moored 
structures  and  vessels,  harbor  resonance,  stability  and  overtopping  of  shore 
protection  structures,  and  surf  beat.  Wave  grouping  is  especially  significant 
in  shallow-water  laboratory  basins  such  as  the  CERC  directional  spectral  wave 
basin.  The  control,  prediction,  measurement,  and  analysis  of  wave  groups  are 
necessary  for  CERC  to  fulfill  its  mission  of  advancing  the  state  of  the  art  in 
coastal  engineering  and  laboratory  physical  modeling. 

52.  Successive  wave  heights  are  dependent  phenomena.  Thus,  the  Kimura 
model  is  a  better  predictor  of  run  lengths  than  the  Coda  model.  The  computer 
program  KIMUR5  gives  excellent  agreement  with  van  Vledder's  values.  Addi¬ 
tional  development  and  testing  with  simulated  and  measured  data  might  lead  to 
better  agreement. 

53.  Future  enhancements  might  include  converting  the  main  program  to  a 
subroutine  so  that  it  can  be  called  from  a  SIWEH  analysis  and/or  spectral 
analysis  program.  Presently  the  correlation  coefficient  is  an  input  param¬ 
eter.  An  option  could  be  to  calculate  it  directly,  from  the  wave  height 
time-history  using  an  autocorrelation  procedure,  or  spectrally,  using  Battjes' 
or  Goda's  method.  Finally,  the  Rayleigh  PDF's  are  calculated  using  the  mean 
wave  height.  An  option  to  allow  the  use  of  other  wave  height  values,  such  as 
median  and  significant  wave  heights,  could  be  included. 

54.  Correlation  among  the  spectral  peakedness  parameter,  the  correla¬ 
tion  coefficient,  and  the  groupiness  factor  obtained  from  a  SIWEH  time-history 
could  be  investigated  further.  Also,  additional  research  into  the  relation¬ 
ship  between  the  groupiness  factor  and  the  degree  of  grouping  in  simulated 
data  and  wind-generated  waves  would  be  beneficial  to  increase  our  understand¬ 
ing  into  wave  grouping  physics. 

55.  An  analogous  program  could  be  developed  for  the  groupiness  statis¬ 
tics  of  successive  wave  periods  which  fall  within  a  certain  period  band. 
Usually,  values  between  0.7  to  1.2  times  the  mean  wave  period  are  most  impor¬ 
tant.  The  development  for  wave  periods  again  assumes  the  time  series  of 
periods  is  mutually  correlated  and  forms  a  Markov  Chain.  Kimura  (1980)  showed 


that  the  Weibull  distribution  would  replace  the  corresponding  one-  and  two- 
dimensional  Rayleigh  distributions.  The  spectral  width  parameter  was  the 
analogous  spectral  parameter  (i.e.  spectral  peakedness  for  heights)  most 
closely  associated  with  the  correlation  coefficient  for  wave  periods. 
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APPENDIX  A:  LISTING  OF  PROGRAM  KIMUR5 


PROGP AH  KIRUR5 

SINGLE  PRECISION*  RANUAL  INPUT*  RIHENSI ONED  FOP.  50  1  ARRAY  SIZE 

PRCGRAH  TO  CALCULATE  GROUP  RUN  LENGTH  STATISTICS  9  ATE  P  ON  KIHURA'S 
HETHOO.  THE  TIHE  SERIES  OF  WAVE  HEIGHTS  ARE  ASSURED  TO  PE  HUTTUALLY 
CORRELATED  AND  FORH  A  HARKOV  CHAIN.  TRANSITION  PR CHAB I L I T  I ES  ARE 
CETPRHTNED  FROR  THE  TWO-DIHENSIONAL  RAYLEIGH  JOINT  t-ROE*PILITY  DENSITY 
FUNCTION  FOR  SUCCESSIVE  WAVES.  FOR  BOTH  A  PUN  OF  UGH  WAVES  ANO  A 
TCTAL  RUN*  THE  PROBABILITIES  OF  RUNS  OF  DIFFERENT  LrNG  THS*  THE  AVERAGE 
RUN  LF.NCTH,  ANC  THE  STANDARD  DEVIATION  CF  THE  PUN  LFNGTH  ARE  CALCULATE 

REAL  Q<E01>»PHR<25)*PTR<2S1 
Rr  A  L  K «  J1 R  « J2H 
PEAL  PPfOltEOlJ 

DATA  R HH1 «NH *r E LH* HH *HC/.6R 1 *97* . 0567 *2 . 1*2 . 1 / 

OPEN  PTSK  FILE  FOR  OUTPUT 

CPENIUf.  TTr2,FILr  =  »K  1HUR  .OU  T  «  «  S  T  ATUS  =*UN  KNOW  N*  ) 

CU'FY  U5EF  FCF  TEST  CASE 
k'R  I  TE  ( f  *  1  C  ) 

FOR  RAT  C  /  *  *  FNTF.p  1  FOR  TEST  CA  St:  »> 

REACT‘D**)  I  A  S  K 

IF  <  I A  SK  .EG.  1)  GO  TO  20 

INPUT  P.PARETFPS 

CALL  INPUT<FHH1,NH*L£LH*HH*HCI 
CALCULATE  KAPPA  CORPELATION  PARAHETER 

CONTINUE 

CALL  ft FPA(PHHI,K» 

CALCULATP  RAYLPIGH  PROBABILITY  DENSITY  FUNCTION  G(M> 

CALL  RP IF (NH,CFIP,HR,C) 

CALCULATE  RAYLEIGH  JOINT  F'RC3  ABILITY  DENSITY  FUNCTION  F<H1»H2> 

CALL  R JPCF <Nh*CFLH*HH*K*P> 

CALCULATE  "'AVLfclCH  TRANSIT  1 1.  N  PR  Gfc  A  B IL  !  TIE  C  Ml  f  K.2 
CALL  PTF ANS,rELr*HC*Q.F»Fl!.P2?) 

CALCULATE  HIGH  RUN  STATISTICS 
CALL  HP  UN  <?r»*F  ‘•R*J1R*SI5J1  ) 

CALCULATE  TCTAL  *UN  STATISTICS 
CALL  TFU*.(Pll*P?2*PTR.J2H,SIGjr> 

OUTPUT  RESULT  5.  Tf  T  E  ®  HI  A'  A  L  .  R  LINE  FR  I  f  TF  R 

CALL  CUTF  T«  =  Hhl  ,NH«[:ELH*HR»rC»K«G*P*Pll  *Hr2  *PHH*  vl  ".S,  It -wl  • 

PTF  ,u2R*S ICj;  ) 


STf'P 

END 

SUP  POUT  I NE  INPL  T  C  P  HH1  .NH*CELH»hM»HC> 

QUERIES  USER  FOP  INPUT  PARAMETERS  FOR  KIMURA'S  CRGUF  LENGTH  PROGRAM 
WRITEC6.1 01 

FORMAT  <  /  « •  ENTER  VALUES  FCP:  *, 

A  /«•  RHh  1  rORPELATTON  COFFFICI  '"NT  •» 

A  /♦•  NH  TOTAL  »  CF  WAVE  HEIGHT  MEASUREMENTS  »« 

S  /»•  DEL H  WAVE  HEICHT  INCREMENT 

A  /«»  UPPER  WAVE  FC-KT  INTEGRATION  LIMIT  =  NH  *  DEL  H*» 

A  /»•  HM  MEAN  «AVE  HEIGHT  », 

A  /t*  HC  THPEEHOLI'  WAVE  HEIGHT  ») 

REACT?. *>  RHH1  ,r  H,OELH.HM,hi; 

RETURN 

ENC 

SU. ROUTINE  KAPF'.(PHI,K) 

CALCULATES  CORRELATION  PARAMETER  K  TKAPf  « =G  *PHC  )  GIVEN  C0RrrL#T  ION 
COEFFICIENT  phhi  <  ALSO  Kf.OLS  AS  GAMP  4(H))  U c  I N  C-  SERIES  APPROXIMATION 
METHOD  CF  EATTUES  FOR  THE  ELLIPTIC  INTEGRALS  CF  THf  1ST  &  M.D1  KIND. 

RE  f L  K  *  f I 

INITIALIZE  PARAMETERS 

PI  =  u.  .  A  TAN <  1.0) 

C  r  (If.  -  A.  *  PI)  /  FI 

)  ap r  a  r. 

R  =  C  »  SHH1 
P2  r  P  *  R 

R  3  =  R  *  RS 

K?  r  F,  -  PG/lf.  -  '-?/ 12  ?■ . 

K  =  SCPT(M-) 

RETUF.M 

ENG 

SU‘  FOLTINE  RPCF  'NH.r f LP.HM.r  > 

CALCULATES  1 -D  =AVLEICP  FJf  <  A  E  I  L  I  T  Y  DENSITY  PJNCTIC'  COl). 

IN'UT  VARIAELf 

N M  -  1  CF  !  ME -'VAIS  Or  r  i.  H  -  V  i,.’vr  PM-fT  V  PI-  U  r] 

celf  -  l:lta  sc-.;e‘t  f  dl1'  'ugcls'iv:  *i  w:.v  "ui>i> 

HM  =  MEAf,  WCVE  HEIGHT 

OUTEUT  VAR  TAF  L‘  < 

S  =  l-r  P  *  >  L  E  1 1-  H  FR  0r  Ai  IL  I  TY  MNSI1Y  fiUTICS  .(hi) 

REAL  COM) 


I*  I  T I  A  L  ’  2  F  rCN  '  T  A  N  T  r. 


n  o  o 


P  1 02  =  2.  *  A  T  A  N  ( 1*0) 

Plot  =  f I C?  /  2. 

HM2  =  HM  *  KM 
01  =  P 1 02  /  KM2 
0?  =  P 104  /  HM 2 
MHP1  =  PH  ♦  1 

RAYLEIGH  POP 

DO  I=1«NHP1 

HI  =  <I-1>  *  CELK 
H12  =  HI  *  HI 

CtT)  =  ( Q 1  *  HI)  «  EXP1-  G2  *  K 1 2 ) 
END  DO 


r 

r 

c 

r 

i 

r 

C 

r 

r 

c 

r* 

c 

c 


c 

r 

r 


c 


RETURN 

END 

SUt  POUT INL  F  JPf CCNH*DELH*HM«K*F) 

CALCULATED  2 -D  ‘'AYLEICH  JOINT  5  ROFApILI TV  DtNSITY  FLNCTION  -(K1«H2> 
IN«"L'T  VARIABLES 

NH  =  TOTAL  *  WAVE  HEIGHT  INTERVALS  C'F  TUMMY  VARIAELE  HI  f  H 2 

DEL  H  =  DELTA  INCREMENT  EE  TW'E  r.  N  SUCCESSIVE  WAVE  HEIGHTS 
HM  =  MEAN  UfVE  HEIGHT 

K  =  CORREl »  T  I C  N  KA^A^ETER  KAPPA 

OUTPUT  VARIABLES 

P  =  °  A  YLE I  0 H  JOINT  FEE  P(H1*H2> 

POUHLE  FHETISIGA  iMblfl) 

REAL  e «  E01 tECl  ) 

RE  !  L  K 

INITIALISE  PARAMETERS 
P I  0  a  =  ATAN(l.t) 

P 1 0 2  =  2.  *  PTC* 

p  7  =  a.  »  r-Tp4 

P I  20 A  =  PI  •  PICA 

HM  2  =  HM  •  HM 

KM  4  =  H  *  2  ♦  hM  2 

CH2  =  1 .2  -  K  «  K 

A  1  =  F 12  0*  /  (HM*  .  C  K  2  1 

HI  r  tin.  /  <hm?  •  C  K  2  ) 

XI  -  (P IC2  *  K)  f  ( HM2  •  CM) 

NHP1  =  f,H  ♦  1 

RAYLEIGH  JOINT  FROb A  b  IL1TY  .'■ENDITY  FUNCTION,  P(H1,F2>  =  A  «  F.  *  C 
DC  I  =  1  ♦  H  p  1 

P UMM  Y  HI  WAV?  HEIDHT  VARIABLE 

*-i  =  <i-i)  *  :elh 

H12  =  HI  •  HI 
DO  J=1«NHPJ 

CUMMY  H  2  UAVl  hf  I  (.  H  T  V  *  F  I  A  T  IE 
H2  =  (J-l)  •  f ELH 
H  ? 2  r  nD  • 

H1H2  -  HI  »  H2 
A  F  A  C  T  C  r 


A4 


A  =  *1  *  H1H2 
C  B  r  AC  T L  ? 

B  =  EXC<-  tl  •  <  H 1  2  ♦  HS2>  ) 

C  X  A  F  GURENT  FOP.  ROCIFIEL  Pf.SSEL  FUNCTION  OF  ZEP  C  ORDER 

X  =  XI  •  H1E2 
OX  r  D  P  L  f  (  X  ) 

C  ROCIFIEC  BE'SEL  FUNCTION  OF  ZERO  OR  f  EP 

CALL  eESI<r«OX»PI*l) 

C  FACTOR 

c  =  exf  <x>  *  mi) 

?  FORMAT  C2I5.BF10.3) 

C  RAYLEICH  JOINT  PCF 

P  (  I  *  J  >  =  A  *  F.  *  C 

C  WRIT1  f£»E>  T.J,«.b,X.tin>*C*P<I.J» 

f*d  rc 

EM  OC 

r* 

w 

p  f  t  up  *■ 

CM 

C 

su'-fout t r.r  fesi  »\»v.r i.logt 
l.  VERSION  l.f  -  r  f  crc.  ISP* 


S  1.1 r  F  0  U  T  T  N  E 


CALCULATES  THE  VALLE  CF  THE  BESSEL  FUNCTION  I«  JEX F  (- X  '  *  1  ( K  ,  X » . 


SU.  FCUTIM  CESCF  iPTICN 


THIS  SU  PCUTIHE  CALCULATES  THE  VALUE  OF  THr  FUNCTION 
TFXPI-X  )*I(K,X ) *  WHf  °E  I  IS  THE  FESSEL  FUNCTION  Cr  NON  NEC 
INTEGRA)  (' 6  L  F  R  i‘  =  0.’"t  Ar.C  f  ON  *ICATIVE  REAL  ARGU*'E  *  T  X.  Al 
CALCULATIONS  A  ■-  F  IN  COL'fLl  PRECISION. 


E  e  A  1 1  VC 
Al  L 


LEAST  N  ♦  1 


SUBROUTINES  AND  FUNCTIONS  CALLED: 


NAME  uESCFIPTION 


ERF  OR  CCDES: 


AN  APPROPRIATE  MESSAGE  I'  WRITTEN  TO  THE  LOG  WFtEN  KO  OF  N< 


SUBROUTINE  C  R  E  A  7  I C  N  CATE  ANL  AUTHOR  : 
VF°MGN  l.E  -  Cl  S  E  c  T  1 S70 

-  NFC  CCMFUTATICN  .CENTRE 


SUBROUTINE  PCD IFICAT  ICNSt 


VEF'ICf 
1  .B 

CEE  CF  ifticn: 

1  .(• 


C  '  TE 

r;  Si  p  t  is  to 

r.Pic  i';»i 

O  r  DEC.  1  S  A  4 


AUTHOR /FI D  M 

NRC  COMPUTATION  CENTRE 


GOF  PON  FOI'C-EmE  .'CMP 


crscRjf  tjcn:  Thf  SIT- F  C l T  I  fr  W^r  UPDATE1  to  reflect  current 
GSCA?-  f'C.GRAM  STAf  D  A  F.  D  S  • 


IMPLICIT  rCUbLt  PRECISION  (t-H.C-Z) 

PE-’L  T.<C-FCCT 
CIMENS  I  CN  rid) 

OAT  A  P‘'TN/C.53S760E/ABt'TAK-37/»AMAX/l.ps2e.73  42  779  7r*lF/ 


AMAX»0MIfjz0.c9ScSS999P999P,p-19 
CHECK  THE  A  F  CU '  t  NTS 


IE!*.GE.d  GO  TO  l 
*RITC<LCG,  1  CO)  N«X 

ICO  FORMAT  (/I  0>  t  »Pr5  I  NEGATIVF  ARCLMFN'T.  ORDER 
A  C’T.lt//) 

C I  M )  =  n 
RETURN 
1  N1=N»1 

IFI X.GT.l .p-F >  r  r  to  f 
IF  *  X  »GE  .O.CG)  SO  TC  3 
WRITE(LCC-.l'C)  ,y 

co  r  i  =  i * m 

r  e  i  c  i )  =  t 

RFTl'RN 


'.111** 


ARGUMENT 


A6 


jv  in  jn  •f'0  . "»  > 


'  VJf  A- 


o  o 


C  I  (  K  .  X  )  =  IX/2)»*K/(K  F ACTOF  I  AL  '*  FOR  X  <  ■  l.HP-P 

OEXF<-X)  =  l.PC-X  FOR  X  <=  1.0T-8 

3  ei <i>=i.po-x 

IFIN.EG.P)  RETURN 
ISW=0 

IFIX. EG  .0.001  I GU=  1 

a=o.co 

IF!  ISW . EG. 0  1  A=2.PO«OKIN/X 
CO  5  1=1, N 
EI C  1*1 JrC.Cn 
IF  (ISW.EQ.l  )  GC  TO  S 
IF  (PI  II  )  .GT.AM)  GC  TO  A 
ISW  =  1 
GO  TO  r 

a  E I « !♦!> =p 1 1  n*x /<r*i > 

E  CONTINUE 
RETURN 

r  CALCULATE  STARTING  POINT  FCE  BACKWARD  RECURRENCE 

C 

t  T  =  X 

SOROOT=SGRT (T) 

N2  =  fl 

IF ( T.GT .2P.C  1  GC  TO  7 

M=1  0.e?  *T-R.l  WPS  I  T-C.23  1-1. 2  3*  AF  SIT-1  .7?) -C.A7*A  Ri<  T-F.7E  >♦  1A  .11 

CC  TO  R 

7  ErSPROC  T*  l-P.F'  3C1E  F-2*#E?IT-lc.Of  A39>,''.feP  3R97E-2  • 
t  AFS(T-aa.R232f  )*C  .13«7FAA£-2«.*.FS<T-131  .'»50A»*0.r33A472E-A« 

R  A-.  SIT -IT  A7.E.GFC  >  ♦  0 .32R3 1£  IE -6«  AP  S I  T-230  Gl.  1  BW  .98  1 61  ♦  1  .0 
E  IFIH.GF.Nl  CO  TO  IE 
M  =  N 

IP  #T.GT  ,1A  .c  )  r<  TC  <5 

vi2A71f  .33«7-2Air,7.S  W:'EIT-C.10  5E-31-r  3G.0E  H-O  .51  Af  -2)  - 
A  5C.lt *Af  S (T-C.CC  11-1 3.99* SpS!  T-0 .251-3 .V7* f  ESC  T-0.55  1-1.77- 
«■  AES(T-C.cf •  CESIT-8.R°1'*A1.5 
GC  TO  ir 

«■  "FAXsSC.r  C  CT  *<-C  .«  A7°3Ef  f-l»  IT-1A.AC017)  ♦O.EL3AA?  IE  -1* 
l  tr  <•(  T-AA  ,7Rlr'r  J ♦  C  .F  t71P53r-2*«f  S  IT -167 .3 A3  1*0. 771  7PASE  -3* 

L  AT  S(T-f  3r.27A)*p.r6851c3t-r«A:EIT-351E'.A3  >*0.12A ACA1L 
C  AE SIT-2fFc GP.C )*22.A775)a1.0 
1C  IF  fKMAX.GT.N)  C,  f  TO  12 
DO  11  I=FRAX,N 
11  r,I«IMl=0.''0 
R  =  f  VAX-! 

* :=rrax 
ir 
c 

C.  CALCULATE  T  F  E  F,'TIC  II  P,X  »  /  I  I  H- 1  *  X  > 

r 

L=2.0»SCFCOT*f.f 

RAIIO=r.LC 

TO  13  1=1. L 

FLOT=2* IFaL-I 1 

RATIO:  »/<Fi CTvyoRATIO) 

3  3  COM  I  PUL 

(  CCFFUTE  F  I P  1  ♦  F  I  v - 1  )«•••«  F <0)»  AND  ALPHA 

r 

i  =  1  .r  -l  s 


A7 


ooorr>Torir> 


XX=2.D0/X 
FP*2  -A 

FH1 =  A/P  AT  10 
IF(H.GT.N)  GC  TC  14 
BI<PM)=FH2 
BI<*>=F*1 
GO  TO  15 

14  IF(Hl.EQ.N)  f  I  IMJ-FMl 

15  ALPHA=FM1‘FH2 
00  16  I  =  1  •  F*  1 

FNO=MI‘ XX*FH1*FH2 
IF<«I.LE.N1>  ei(f*I)=FFO 
ALFHA=ALpFA‘FHO 
FM?=FH1 

16  FH1=FM0 

ALPHA  =  2.DOALPHA-FH0 
C 

C  CALCULATE  THE  VALUES  iF  PC  *M -X  »  *  1  IK*  X  >  *K  =  0  »N  • 

r 

IFCALPHA.LE.AMAX1  GO  TC  1* 

A=LHIN* / LFhA 

17  IF »pI CN?) .GT.A  )  GO  TO  IP 
BI«N2)=P.CC 

N2=N?-1 
CO  TO  17 

16  ALPHA=1 .C C /ALPHA 
DC  19  I  - 1  * N2 
19  61 <  I)=B  I<  I) ‘ALPHA 
RETURN 
END 

SUSeCUT  IKt  FTP  »t  H.DCL*"  »HC  «r  «P«rn  »s02> 

CALCULATES  PAYLEIGH  TRANSITION  PR OP  AC  I L  I  T I f 6  PU  A  P2P. 

I N r  UT  VARIABLES 

NH  :  »  CF  INTERVALS  CF  rUPPT  WAVE  HEIGHT  VAFUfLE  Hi  t  12 

CCLH  r  ii  LTA  I NC  FL  if.  N  T  rLTWEEN  SUCCESSIVE  HI  t  HE  WAVE  HEIGHTS 

HC  =  CUTOFF  CR  THFESHCLC  WAVE  HEIGHT*  ALSO  H» 

G  =  RAYLEIGH  1-0  PROBABILITY  TENSITY  FUNCTION  r.CHl> 

P  =  PAYLEIGH  2-0  JOINT  PACE  ABILITY  DENSITY  FUNCTION  PIF1,H2> 

C 

C  OUTPUT  VARIABLE* 

r  P  1  3  r  TRANSITION  PPOE  4[-  IL ITY«  N  r  I  T  H  f  R  HI  NCR  H?  EXCEEDS  THRESHOLD 

C  WAVE  HEIGHT  HC 

C  F 22  =  TRANSITION  PROBABILITY.  60TH  Hi  R  bp  E  X  CE  E  E  THF  ESELLE  WAVE 

C  WAVE  HEIGHT  HC 

C 

REAL  C«rCl> 

REAL  pf*01«5CI> 

C 

C  Pll  TRANSITION  E ROE  ABILITY 

r 

NL  =  1 

NU  =  Hr.  /  DELH 

CALL  RTH  <NL«M*DELH.fi.»  »  f  1 1  ) 

r 

C  P2?  TRANSITION  PPCtABILITY 

C 


A8 


NL  =  HC  /CELH 
NU  =  NH  ♦  1 

CALL  RTFI <NL.NU.OELH.C.F.P?r> 

RETURN 

END 

SUBROUTINE  RTP I (NL .NU .CFLH .G »P  .P90B > 

INTEGRATES  1-C  &  2 -D  9AYLFICH  FPOBABILITY  DENSITY  FUNCTIONS 
Q  (  H 1  )  &  P ( HI »H2  > 

input  variables 

NL  =  LOWER  INTEGRATION  UNIT  ARRAY  ELEPFNT 

NU  =  UPPER  INTEGRATION  LIMIT  ARRAY  ELEMENT 

OELH  r  DELTA  INCREMENT  PETWEEN  SUCCESSIVE  HI  OR  h?  WAVE  HEIGHTS 

P  :  RAYLEIGH  2-0  JOINT  PRCEABILITY  DENSITY  FUNCTION  PCH1.H2) 

Q  =  RAYLEICH  I-C  PROBABILITY  DENSITY  FUNCTION  CCM1» 

OUTPUT  VARIABLES 

PROF  =  TRANSITION  PROBABILITY.  EITHER  Ell  OP  P22 
REAL  P(501.F01> 

REAL  Q  (  E  0  1  ) »PH 1 (501>»PH2(f01>.2<501> 

2-C  RAYLEIGH  JOINT  f  R  0  B  A  6  I L  ETY  DENSITY  EUNCTICN  INTEGRAL 
IN  NUMERATOR 

M  =  0 

SHIFT  PDF  TO  DUMMY  ARRAY  PH  2 
CC  I=NL»NU 
K  =  r 

CO  U=NL.NU 
N  =  K  *  1 
PH2 (N)  :  P(  I.J) 

END  DO 

IE  TEG  FATE  IN  P2  DIRECTION  &  CREATE  NEW  DUMMY  ARRAY  PHI 
NC IM  r  \ 

CALL  GSF (LELh.PHD.Z.NCIM) 

•>  =  F  ♦  ’ 

F  HI ( M  >  =  7 ( NO  I M ) 

FNC  00 

INTFGRATE  IN  HI  DIRECTION  USINC  DUMMY  ARRAY  PHI 
NO  I M  =  v 

C*LL  QSF(rELh*PH1.2.NriF> 

SUMN  =  7(N:iM> 

1-0  PAYLEIGH  PROBABILITY  CENSITY  FUNCTION  INTEGRAL  IN  OENCMNATCR 

SHIFT  EPF  TO  DUMMY  ARRAY  PHI 
N  =  0 

00  I  =  M  .  NL' 

N  =  N  ♦  1 
P  HI <  N  )  =  G  < I  ) 

ENC  DO 

INTEGRATE  IN  HI  C I  RFC  T  I  f N  USING  DUMMY  ARRAY  PHI 
NC I M  :  N 

CALL  CSF  CDELH»rHl»Z*NPIM) 

SUMC  =•  7  C N  D  IM  > 


o  o  r>  o  o  r*>  r>  o  o  r>  o  ^  n  o  o  o  o  o  o  r>  o  o  n  o  o  o  o  <■> 


! 


TRANSITION  RROc  A6 IL  !  T Y 

PRCB  =  SUNN  /  SI'*»C 

RETURN 

END 

SUEROUTINE  OSF I R« Y « 2 »ND  IP ) 


SUBROUTINE  «SF 
'’UR°OSE 

TG  COMPUTE  THE  VECTOR  OF  INTEGRAL  VALUES  FOR  A  GIVEN 
FCUILISTaNT  TABLE  OF  FUNCTION  VALUES. 

USAGE 

CALL  OSF  <H,Y*Z*NriM> 

tescwipticn  r r  parameters 

►  -  THE  I NC  SFMr NT  0^  Afi  f-UMEN  T  V  i  LUES . 

Y  -  THE  1NFLT  VECTOR  CF  FUNCTION  VALUES. 

7  -  THE  CE  SULT I NC  VECTOR  OF  INTEGRAL  VALLES.  2  NAY  EE 

TTENTICAL  WITH  Y. 

NTIN  -  THE  r.INENSION  OF  VECTORS  Y  AND  2. 

PE  HA  0  K  S 

NC  ACTION  IN  CASE  NCIH  Lt'S  THAN  3. 

SUEROUTINES  ANC  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

.w  E  T  u '  V 

FfGINMNC  WITH  ?<l)=r,  EVALUATION  OF  VCCTEF,  t  IS  ICNf  U  Y 
“FANS  OF  SIMPSONS.  R  UL '  T  C  C  EThER  WITh  RFWTONS  */f  RULE  OF  ^ 
CCPIINATIfN  CF  THE Sr  TWO  PULES.  TPUNC-’TION  ERROR  IS  CT 
ORDER  H*  • 1  II. E.  FOURTH  ORDER  METHOD.  ONLY  IN  CaSL  NriN=3 
TcU»CATIiV  E c. 0  G  R  CF  7'2)  IS  OF  OR:rP  H««*. 
cff  RfFE  ••ENCE  .  '"EE 

II)  F.E.I  ILCEEKANr'»  I  NTPC  DUCT  ION  10  NLMERICAL  ANALYSIS* 
"CCRAW-HILL*  NEW  YC»K/TOPOfUO /LONEON.  19E.fe*  T R. 71-76. 

C  t2)  c • 2UR “Ur  RL  *  PRA*TISCK  MUHEMATIK  FUER  I R  Ot  NI  EURE  UNO 

:  RBYSlRER.  SPPINGtR.  trC  RL  I N/GC  F  TT  I  NGf  \  /HE  I  f.  E  LF>  E  P  G  .  1EE3* 


C 


f 


r 


D  I  FENS  ICN  Y(P.7(I) 

HT=.?3J.*?3?«R 
TF  I  Mr  I 

NdP  IS  GREATER  TRAN  E.  PREPARATIONS  OF  INTI  GRiAT  ION  LOCR 
SU*  1  =  V  I  ?  )  ♦  Y  I  2  ) 

SUM’rSUMI  ♦SLA*  I 

SUN?  =hT*<  Y<  1  I«.LA  1*Y  <  *  >  ) 

AL'VIrYIA>*NC4) 

AUX1=AUX1*AUX1 

AUY1=SU*1*R  T»I YI 3)*SU»1*YI' J  > 

AUX',=HT*IYn)»T.R7b*iYI2>  ♦»(*>»  ♦2.G?E*IYl3J*Y(A>)^tif)) 
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SUf*?  =  V<5J»Y'?> 

$Uf*2=SUN2*SUW2 

SU*2=AUX2-HT*<  V  (M*SUN2*Y  (O  ) 

2<1)=0. 

AUX=Y(3  )*Y(3) 

AUX=AUX*AUX 

2<2>  =  SU*2-HT*( Y(2>»AUX*Y<4>  > 

2<?>=SU*1 
2  (4  >=SUP2 
IFfNDIN-fc>5»5*2 
C 

C  INTEGRATION.  LOOP 

2  DO  4  1  =  7*  NO  IN*  2 
SUM1=AUX1 
SUP2=AUX2 

AUX 1  =  Y  < !-l)*y(I-t) 

AUX1=AUX1*AUX1 

AUX  1  =  SUW1*HT*(Y (!-2>*AUXl»Y<I>> 

2(  I-2J  =  CUX*1 
I F  {  I-Nnif‘)7»fc»f 

3  AUX2=Y<  I  )  ♦ Y  M ) 

A  U  X  2  =  A  U  X  2  ♦  *  U  X  2 

AUX?  =  SUf?*hT»( Y ( I-l)>AUX2*Yf  I  ♦  I  >  > 

a  2  ( i  -l  >  =  fu**z 
5  7(ND1H-1»=AUX1 
2<NCIN>=ALX? 

Bf  TURN' 

F  Z(NOI"-l)  =  SU*2 
2<N?,INI  =  AUX1 
BF  TURN 
C 

0  ENT  OF  JMtGPATIOA  LOOP 

C 

7  IF(NPI*-3U2*11*F 
C 

C  NEIP  I?  ECUAL  TO  4  OF  * 

e  SUM2=1»125*HT*( x  (  1 >*Y<2)4Y<?>4YC2)«Y(2}*Y(3>*Y(3>*Y(A  >) 
SlVl=Y(D*Yf2I 
SUM  =  EUM  ♦SUf'l 
EUF  1  =  1  T«(Y(1>*MIP1*Y(,»> 

Z<1 >=n. 

AUX !=V C’)yY<3> 

AUX  l  =  Ai:  X  1  ♦  A IIX 1 

2<2  >  =SLf  2-*-  T«  (  Y  <  2  >»£UX  1*Y  «  4  U 

IMF.r  I“-‘.  M0**»tV 

c.  Al’Yl:Y(A)«VIA) 

AUX 1  =  AU  XI *AUX1 

2(5>  =  SUP1'»uT*(V(3>*AL,X1yY(->) 

10  2(3>=S UM 
7f4)rru'  2 
0  E  T  UF  N 

C 

C  NO  I  *  TS>  ECUAL  U.  7 

11  SUNl=HT*(l.re*Y<l)*Y«7>*Y<?>-.2E»Y(;n 
<;U»'?  =  Y«?)*Y(?) 

SU"?  =  SU*’2»SU*? 

2(1  >  =  0. 

Z(2):SI'M 

12  P  F  T  UP  N 


All 


SU '  ROUT  ?N»  HRUMF22*FHR,JlM«SIC-ul) 


C ALCUl A  TES  HIGH  RUM  GROUP  STATISTICS  OF  PROBABILITY  FOR  DIFFERENT 
RUN  LENGTHS*  ME AN  RUN  LENGTH*  t  STANDART  DEVIATION  CP  RUN  LENGTH. 

INPUT  V  AP  I  ABLE  S 

P02  -  T  F  ANS  !  T  !  CM  PPOFAPUITY  FCR  SIHLLT  ANEOLS  E  >  CE  EG  ANCE  LF 

THREShCLC  FT  fcOT*  HI  t  H2  LAVE  HEIGHTS 

OUTPUT  VARIABLES 

FHR  -  FROEAEILITV  OF  RUN  LENGTH  HAVING  LENGTH  OF  Jl«  F  I «.  1 ) 

JIM  r  M  f  A  N  M'N  LENGTH 

S  I  r-  J !  =  STANC’FPC  Lc  V  IATICN  OF  ‘UN  LENGTH 

REAL  Jlf> 

REAL  PhF  f  TS  ) 

P°C r  AE  I  1  I T  Y  Of  MIN  OR  LfNGTH  Jl.  PHKCJl) 

IF  <  r  22  .EC.  O.T  °22=l.r-0F 
Cl • 22  -  1."  -  tit 
UO  J1=1.?E 

FfPEJI)  r  P2;**(tl-1)  *  r 1 p  2  2 


EM 

CC 

ME  AN  Rt; 

H  L:l.fTH 

JIM  =  1 

.  /  r.i»rf 

ST  *' 

.I  A  r 

"  r  r  V  »  A  T  Tfr  r  f  P(iN 

' Iv  J1  = 

SCFTAF;  ?)  /  C 1 P?2 

F  F  T  UP  N 
f  >:r 


SU;  F  C  U  T  IN'  TRU".  .'PI  1  *F2  2«FTF  *jrf*.5lGJ2) 

CALCULATES  TCTAL  pin  CRCUF  STATISTICS  F  PROBABILITY  FOR  CIEFCRENT 
FUN  LTNCTHS,  mean  hUN  LENGTH*  S  STANDARC  DEVIATION  (F  F  LN  LENGTH. 

IDF  L  T  P  I  '  •  l *  E 

c  !  1  -  TRANSITION  PFC^A'IIITY  FDR  NFITHEF  HI  f  OF  M  F  >  C*. !.  t  1  * '• 

TF  R  E  S  *-DLL  LAVE  HEIGHT  1  C 

(  2  ;  r  T  =  &NS  T  T  I  0  F  PPOEAEIIITY  FCR  SIMULTANEOUS  F  WEE  F  ANCE  OF 
T  t-F  E  SHOLC  HC  FY  BOTH  Hi  l  H2  LAVE  FEIGFTS 

OUTPUT  V  t  F  I  1 1-  L  i  S. 

FTC  r  F  c  C  L  A  i  ’LITY  CF  H'  .  LENGTH  HAVIK  L  r  NC  T  h  CF  J2»  FIj:) 

J2  H  =  MEAN  'IN  L'  NOTH 

SICJ2  =  S.T'NCM-O  DEVIATION  OF  FUN  LF  NC  TH 

FFfi  j;? 

real  r T 1 ( r 1 ) 


f  F‘Ot  ABILITY  CF  ’UN  0  F  LENGTH  J'.  rTFfjr> 

If  *  r  1  1  .E  r  .  0. )  M  lrl . E-C A 


A12 


:^v  v.v.’a 


sC 


IF*r?2  .EG.  C.)  F2r=l.r-o; 

C!  =  1.0  -  rn 
C2  =  1  .0  -  F22 
C*  =  Cl  •  C 2  /  < P 1 1  -  P2C) 

CO  J2=l»25 

PTR1J21  =  C3  •  lPll**»J?-i>  -  P2?»MJ"-1  )  ) 

END  00 

PEAN  RUN  LrNGTF 

J?P  :  1,  /  fl  «  1.  /  r? 

STANDARD  DEVIATION  CF  ^ OK  LEN6TM 

C 1 2  =  Cl  •  Cl 
C22  =  C2  •  C2 

SICJ2  =  SCRT<P22/C:2  ♦  I11/M2) 

RETURN 

END 

SU  BOUTIN?  CUTFT<t‘-Fl»FMtffLM»U*'»FC»K»Ci«r’*Fil»P21*F  R  .JlR» 

S  K  J1  .(  TR  .  J2F  .5  I  C  J  2  > 

OUTPUTS  RESULTS  FFC*  «T»URA‘S  ‘LGDFITHP  FC«  DAlCULAlINr~  UUF 
LEFGTH  STATISTICS  ‘Y  ASSlNi.C  T-tT  LAWS  t-EKMTG  A°  r  CORRELATES. 

RE*  L  P  (rPl.c01 ) 

REAL  G * * Cl  )  .PHK <2* ) »F In < ?’  ) 

REAL  K.U1F»J2« 

OE SCR  IF  TI  VC  TITLE 
;  T»  C  «1C) 

FORMAT*  I*]  ,/»TTF»,rt  SUITE  c  C  *  1  F  C  C  F  A  ►  K  ’  v  u  ,  '  •  » 

✓  *TaL  »•(  RCL'F  Pi*  L  C  A  C,  T  U  STAT  I  ST  I r  S • > 

T»R  UT  V  A F  14«=Lr 

MU  *  '.h  •  !FLU 

r  m  >  :  =  r  m  ♦  i 

WR  I  TF<  ?*?r  >  FHt-1  *NF  *t  !  L  F*MI  .HR«MC 
FCRVAT<//.*  •*  •  .  •  I  Nf  L  T  VAR  lffiLEo«**«»*  « 

//.’  rORRELATIC'  SO!  CF  K  irMTt  SfM  -  •*T*«r»Flr.R» 

/«*  UT'L  ■  C  F  ;AV'  M  f  T  >.  n  T  F  i  A  0  I  M  *•  :  K  T  S  «  \r  =»»74*»,ir« 

/.•  VAVR  *-!  It'-T  If.  LEL*-  *  •  .  T*r  *F  1C  .<•• 

/,»  Ut  *■  E  F  ,  A  V  E  f-CMT  I  N  T 1  ORATION  Li*  I  T  -  •  .  T  <*  ‘  t  F  1  r  .  <•  « 

/.•  F  E  a  f  WVf  MFKfl.  *•»  =*,TAC 

/«•  THF S  SMGLO  WAVE  ' l !'  1  T»  MC  =  •  i TN -  . R 1 C  .  A  ) 

SI  TCLT  V t R  I  A t  L  ;  ' 

•  c  l  tf • - .  :c  >  f . -  u  :  ? 

FOA‘ATI//.»  .....rLTFLT  VAR  1  At LL3*«  ••••  . 

//,•  CC<-  R  EL  /  T  Iff  t- it.  •.PETER  ,  K 

/,•  f«C!  AMI  ITr  '  I  I  Th.fR  Hi  FOR  H?EvrL-fr««  ‘11  =**T<,‘.U0.« 
/«•  RFC'  '■ 1  1 1  ITT  ,'T  f-  1  >.  H'  SXCi'r.S.  F  2  2  =•  »T«.E  *F1  f.A  T 
VR  I Tr  t  -  ,»>  ) 

FM  “'T  ( /,  •  r  A  Y  L  i  I  T- »  1-.,  F'P  V-*L  FS  fKE:*> 

bR  I  Tr  I  '  r  *  H  I) 

F  C  R  **A  T  (  1  r  I  /.■>) 


W«  ITE<2  ,*5> 

FGFFAT</,»  1  A  VL  f  1 0  F  ?-C  f-  CF  VALLES  ARE:*) 

WRITE  12  ,4  0  »  <CP<  I,U>,  I  =  1»T:HP1  >  ,«  =  1  ,f.MPl) 

WR  I  TE  <  2 ,5  C  ) 

FORMAT (✓/*•  *,***FISF  WAVE  RUN  fRCUP  RESULT $••••**  , 

4  /,•  PR 0£ ABILITIES,  THR:») 

WRITE<2»40>  CPHMI  >,1  =  1,25) 

WRITE  <2, 60)  J1H,SI0J1 

FOR  FAT  </«•  *E  A  N  GF  CUP  LENOTF,  Jl*  =  •  ,T  4  *  ,F  1  n.  *  , 

4  /,*  STD  CEV  GROUP  LF'GTh,  SIGdl  =•, T45 ,F 1 C. 4  ) 

wP I TE  f  2  *7  0  ) 

FCR«AT(//,*  ,,.,*TCTAL  WAV'"  RUN  GROUP  R  lS  UL  TS  **•**•  , 
t  /,«  PROBABILITIES,  rTP;») 

WRITEI2,<*r>  *PTMI>, 1  =  1,25) 

WPTTEtS.RC)  J2N  »  SI  G»J? 

F0C  FAT ( /*•  FEAN  CfOUP  LENGTF,  w2M  =  •  ,T  4  r  ,F  1  C  .4  . 

K  /,♦  STD  TEV  ''-POUF  LEFC-TF,  SI6J2  =  •  ,  T  45  *F  1 0.  <•  > 

RETURN 

END 


APPENDIX  B:  LIST  OF  SYMBOLS  USED  IN  PROGRAM  KIMUR5 


Symbol 

DELH 

H 

HM 

HC 

HU 

HI 

H2 

JIM 

J2M 

K 

NH 

NL 

NU 

P(H1 ,H2) 

PHR(Jl) 

PI 

PROB 

PTR(J2) 

Pll 

P22 

Q(H1) 

RHH1 

SIGJ1 

SIGJ2 

Y 

Z 


_ Description _ 

Delta  wave  height  increment  between  successive  HI  or  H2  wave 
heights,  controls  upper  wave  height  integration  limit,  HU 
=  NH  *  DELH  in  transition  probabilities 

Increment  of  argument  values  (i.e.  X-array)  for  calculating 
integral 

Mean  wave  height 

Cutoff  or  threshold  wave  height 

Upper  wave  height  integration  limit 
HU  =  NH  *  DELH 

Dummy  wave  height  variable 

Dummy  wave  height  variable 

Mean  run  length  for  run  of  high  waves 

Mean  run  length  for  total  run 

Correlation  parameter 

Total  number  of  wave  height  measurements  or  intervals  of  dummy 
wave  height  parameters  HI  or  H2  between  zero  and  upper  wave 
height  HU  in  transition  probabilities 

Lower  integration  limit  array  element 

Upper  integration  limit  array  element 

2-D  Rayleigh  joint  probability  density  function  for  successive 
wave  heights 

Probability  of  run  length  having  length  of  J1  for  run  of  high 
waves 

3.14159... 

Dummy  transition  probability,  either  Pll  or  P22 

Probability  of  run  length  having  length  of  J2  for  total  run 

Transition  probability,  neither  HI  nor  H2  exceeds  threshold  wave 
height  HC 

Transition  probability,  both  HI  and  H2  successive  wave  heights 
exceed  threshold  wave  height  HC 

1-D  Rayleigh  probability  density  function  for  individual  wave 
heights 

Correlation  coefficient 

Standard  deviation  of  run  length  for  run  of  high  waves 
Standard  deviation  of  run  length  for  total  run 
Function  values  (i.e.  Y-array)  to  be  integrated 
Vector  array  of  integrated  values 


B2 


APPENDIX  C:  DESCRIPTION  OF  SUBROUTINES 


Subroutine  BESI 

Description:  Calculates  value  of  Bessel  Function  I,  Dexp(-X)  *  I(K,X) 

where  I  is  Bessel  function  of  non-negative  integral 
order  K=0,  N  and  non-negative  real  argument  X.  All 
calculations  are  in  double  precision.  National 
Research  Council  of  Canada  (NRCC)  Scientific  Library 
Subroutine . 

Calling  Statement:  SUBROUTINE  BESI  (N, X.BI.LOG) 

Arguments: 

N  1*4  Order  of  Bessel  Function  I,  N  =  0  for  zero  order  (used 

here) 

X  R*8  Argument  of  Bessel  Function  I 

BI  R*8  I-D  array  containing  value  of  function  Dexp(-X)  *  I(K,X) 

in  its  K+l  element  for  K=0,N.  Dimension  of  BI  in  call¬ 
ing  program  must  be  at  least  N+I,  equals  1  for  N=0  for 
modified  zeroth  order 

LOG  1*2  Logical  unit  number  of  error  message  when  X<0  or  N<0 

Called  By:  Subroutine  RJPDF 
Calls  To:  None 

Reference:  NRCC  Computation  Center 
Subroutine  HRUN 

Description:  Calculates  high  run  group  statistics  of  probabilities  for 

different  run  lengths,  mean  run  length,  and  standard 
deviation  of  run  length. 

Calling  Statement:  SUBROUTINE  HRUN  (P22 ,PHR, JIM.SIGJ1 ) 

Arguments : 

P22  R*4  Transition  probability,  both  HI  and  H2  successive  wave 

heights  exceed  threshold  wave  height  HC 
PHR  R*4  Probability  of  run  length  having  length  of  Jl,  PI(J) 

JIM  R*4  Mean  run  length  for  a  run  of  high  waves 

SIGJ1  R*4  Standard  deviation  of  run  of  length  Jl  for  run  of  high 

waves 

Called  By:  Program  KIMUR5 

Calls  To:  None 

References*:  Van  Vledder  (1983a, b) 

Kimura  (1980) 

Coda  (9185) 


*  References  in  appendixes  are  cited  in  References  at  the  end  of  the  main 
text . 


Subroutine  INPUT 
Description: 


Queries  user  for  five  input  parameters  for  Kimura's  wave 
group  analysis. 


Calling  Statement:  SUBROUTINE  INPUT(RHH1 ,NH,DELH,HM,HC) 


Arguments: 


RHH1 

R*4 

Correlation  coefficient 

NH 

1*2 

Total  number  of  wave  height  measurements  or 

total  number  of  intervals  of  dummy  wave  height 
parameters  HI  and  H2  between  zero  and  upper  wave  height 
HU  in  transition  probabilities 

DELH 

R*4 

Delta  wave  height  increment,  controls  upper  wave  height 
integration  limit,  HU  =  NH  *  DELH  in  transition 
probabilities 

HM 

R*4 

Mean  wave  height 

HC 

R*4 

Cutoff  or  threshold  wave  height 

Called 

By: 

Program 

KIMUR5 

Calls 

To: 

None 

Reference : 

None 

Subroutine  KAPPA 

Description:  Calculates  correlation  parameter  Kappa  given  correlation 

coefficient  using  series  approximation  method  of  Battjes 
for  Complete  Elliptic  Integrals  of  1st  &  2nd  kind. 

Calling  Statement:  SUBROUTINE  KAPPA  (RHHI.K) 

Arguments : 


RHH1  R*4  Correlation  coefficient 

K  R*4  Correlation  parameter 


Called  By: 
Calls  To: 
Reference : 


Program  KIMUR5 
None 

Van  Vledder  (1983a, b) 


Subroutine  OUTPT 
Description: 


Outputs  results  from  Kimura's  algorithms  for  calculated 
group  run  length  statistics  to  disk  file,  logical 
unit  2,  KIMUR . OUT . 

Calling  Statement:  SUBROUTINE  0UTPT(RHH1 ,NH,DELH,HM,HC,K,Q,P,P1 1 ,P22,PHR, 

J1M.SIGJ1 ,PTR, J2M,SIGJ2) 


Arguments: 


RHH1 

R*4 

Correlation  coefficient 

NH 

1*2 

Total  number  of  wave  height  measurements 

Total  number  of  intervals  of  dummy  wave  height  parameters 
HI  and  H2  between  zero  and  upper  wave  height  HU  in 
transition  probabilities 

DELH 

R*4 

Delta  wave  height  increment  between  successive  HI  and  H2 
dummy  wave  heights,  controls  upper  wave  height  integra¬ 
tion  limit,  HU  =  NH  *  DELH  in  transition  probabilities 

HM 

R*4 

Mean  wave  height 

HC 

R*4 

Cutoff  or  threshold  wave  height 

Q 

R*4 

Rayleigh  one-dimensional  (1-D)  PDF  Q(H1) 

P 

R*4 

Rayleigh  two-dimensional  (2-D)  PDF  P(H1,H2) 

Pll 

R*4 

Transition  probability,  neither  HI  nor  H2  successive  wave 
height  exceeds  threshold  wave  height  HC 

P22 

R*4 

Transition  probability,  both  HI  and  H2  successive  wave 
heights  exceed  threshold  wave  height  HC 

PHR 

R*4 

Probability  of  run  length  having  length  of  Jl,  P1(J) 

JIM 

R*4 

Mean  run  length  for  a  run  of  high  waves 

SIGJ1 

R*4 

Standard  deviation  of  run  of  length  Jl  for  run  of  high 
waves 

PTR 

R*4 

Probability  of  total  run  length  having  length  of  J2,  P2(J) 

J2M 

R*4 

Mean  total  run  length 

SIGJ2 

R*4 

Standard  deviation  of  run  of  length  J2  for  total  run 
length 

Called 

By: 

Program 

KIMUR5 

Calls 

To: 

None 

Reference : 

None 

Subroutine  QSF 
Description: 


Calling  Statement: 
Arguments: 


H 

Y 

Z 

NDIM 


R*4 

R*4 

R*4 

1*2 


Computes  vector  of  integral  values  for  a  given  equidistant 
table  of  function  values.  Computes  integral  of  function 
contained  in  array  Y  of  dimension  NDIM  for  equidistant 
X-array  values  spaced  H  apart  using  combination  of 
Simpson's  and  Newton’s  3/8  Rules. 

SUBROUTINE  QSF(H,Y,Z ,NDIM) 


Called  By: 
Calls  To: 
References : 


Increment  of  argument  values,  DELH  for  x-axis  array 
Input  vector  of  function  values 

Resulting  vector  of  integral  values,  contains  integrated 
area  under  curve  represented  by  array  Y 
Dimension  of  vectors  Y  and  Z 

Subroutine  RTPI 

None 

NRCC  Computation  Center 
Hildebrand  (1956) 

Zurmehl  (1963) 


Subroutine  RJPDF 

Description:  Calculates  2-D,  joint,  or  Bivariate  Rayleigh  Probability 

Density  Function  P(H1,H2)  based  on  Kimura's  theory. 
Uses  modified  Bessel  function  of  zero  order. 


Calling  Statement:  SUBROUTINE  RJPDF(NH,DELH,HM,K,P) 
Arguments: 


NH 

1*2 

Total  number  of  intervals  of  dummy  wave  height  parameters 
HI  and  H2  between  zero  and  upper  wave  height  HU  in 
transition  probabilities 

DELH 

R*4 

Delta  wave  height  increment  between  successive  HI  and  H2 
dummy  wave  heights,  controls  upper  wave  height  integra¬ 
tion  limit,  HU  *  NH  *  DELH  in  transition  probabilities 

HM 

R*4 

Mean  wave  height 

K 

R*4 

Correlation  parameter 

P 

R*4 

2-D  Rayleigh  PDF  P(H1,H2) 

Called  By:  Program  KIMUR5 

Calls  To:  Subroutine  BESI  (NRCC  Scientific  Library  Subroutine) 

References:  Van  Vledder  (1983a, b) 

Kimura  (1980) 

Goda  (1985) 


Subroutine  RTPI 
Description: 


Integrates  I-D  and  2-D  Rayleigh  PDF  Q(H1)  and  P(H1,H2), 
respectively. 

Calling  Statement:  SUBROUTINE  RTPI (NL,NU,DELH,Q,P,PROB) 


Arguments: 


NL 

1*2 

Lower  integration  limit  array  element 

Nl' 

1*2 

Upper  integration  limit  array  element 

DELH 

R*4 

Delta  wave  height  increment  between  successive  HI  or  H2 
dummy  wave  heights,  controls  upper  wave  height  integra 
tion  limit,  HU  =  NH  *  DELH  in  transition  probabilities 

P 

R*4 

2-D  Rayleigh  PDF  P(H1,H2) 

Q 

R*4 

1-D  Rayleigh  PDF  Q(H1) 

PROB 

R*4 

Transition  Probability,  either  Pll  or  P22 

Called 

By : 

Subroutine  RTP 

Calls 

To: 

Subroutine  QSF  (NRCC  Scientific  Library  Subroutine) 

References : 

Van  Vledder  (1983a, b) 

Kimura  (1980) 

Goda  (1985) 

Subroutine  TRUN 

Description: 

Calling  Statement: 

Arguments : 

Pll 

R*4 

P22 

R*4 

PTR 

R*4 

J2M 

R*4 

SIGJ2 

R*4 

Called 

By: 

Progr 

Calls 

To: 

None 

References : 

Van  V 

Calculates  total  run  group  statistics  of  probability  for 
different  run  lengths,  mean  run  length,  and  standard 
deviation  of  run  length. 

SUBROUT INE  TRUN (Pll , P2 2 , PTR , J2M , SIGJ2 ) 


Transition  probability,  neither  HI  nor  H2  successive  wave 
heights  exceed  threshold  wave  height  HC 
Transition  probability,  both  HI  and  H2  successive  wave 
heights  exceed  threshold  wave  height  HC 
Probability  of  total  run  length  having  length  of  J2,  P2(J) 
Mean  total  run  length 

Standard  deviation  of  run  of  length  J2  for  total  run 
length 


Kimura  (1980) 
Goda  (1985) 
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APPENDIX  D 


NOTATION 


E[  ] 
E(  ) 
f 
H 


H. 

l 

*i+l 

H 


1 


H„ 


in 


H 


H 


max 

i 

med 

H 

r 

H 


V  1 

jl 

h 
h 
j  2 

k 

K(  ) 


m 

o 


N 


P 


P(Hj,H2) 

Pn 

po 

pap 

P(j2) 


Expectation  operator 

Complete  elliptic  integral  of  the  second  kind 

Frequency  variable 

Wave  height  variable 

Dummy  wave  height  variable 

Successive  wave  height  variable 

First  of  two  successive  wave  heights  dummy  variable 

Second  of  two  successive  wave  heights  dummy  variable 

Cutoff  or  threshold  wave  height 

Mean  wave  height 

Maximum  wave  height 

Median  wave  height 

RMS  wave  height 

Significant  wave  height 

Cutoff  or  threshold  wave  height 

Modified  Bessel  function  of  zeroth  order 

Run  length  for  run  of  high  waves,  i.e.  1 , 2, 3, ...  1 1 , 12+ 

Run  length  for  total  run,  i.e.  2 ,3,4, . . . 1 1 , 12+ 

Mean  run  length  for  run  of  high  waves 

Mean  run  length  for  total  run 

Lag  of  autocorrelation  function  estimate 

Complete  elliptic  integral  of  the  first  kind 

Zeroth  moment  of  time  series  of  wave  elevations 

Total  number  of  points  in  wave  height  time  series 

Probability  that  wave  height  H  exceeds  threshold  wave  height  or 
Markov  Chain  transition  probability  matrix 

Transition  probabilities — neither  H.  nor  H.  exceeds  threshold 
height 

Transition  probabilities — both  and  exceed  threshold 

height;  simultaneous  exceedance  of  threshold  wave  height  by 
both  Hj  and  waves 

Joint  or  bivariate  Rayleigh  probability  density  function 
Markov  Chain  distribution  after  n-time  transitions 
Initial  Markov  Chain  distribution 
Run  length  probability  for  run  of  high  waves 
Run  length  probability  for  total  run 


D2 


q  Probability  that  wave  height  H  does  not  exceed  threshold  wave 
height 

Qp  Goda's  spectral  peakedness  factor 

q(Hj)  Rayleigh  probability  density  function  for  individual  wave 
heights 

R,,  (1)  Correlation  coefficient  for  successive  wave  heights 

nh 

S(f)  Spectral  estimate  of  the  surface  elevation 

T  Mean  zero-crossing  wave  period 

m 

<  Correlation  parameter,  equals  2p 

p  Correlation  parameter,  equals  k/2 

7i  3.14159.  .  . 

Correlation  coefficient  for  successive  wave  heights 
o  Standard  deviation  of  wave  height  time  series 

o(jj)  Standard  deviation  for  run  of  high  waves 

oCj^)  Standard  deviation  for  total  run 


D3 


